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PREFACE 

The  work  reported  herein  was  performed  by  the  Arnold  Engineering  Development  Center 
(AEDC),  Air  Force  Systems  Command  (AFSQ  at  the  request  of  the  Directorate  of  Technology 
(DOT),  AEDC.  The  DOT  Project  Manager  was  Robert  W.  Smith.  The  Calspan  Project 
Manager  was  William  T.  Strike.  The  results  were  obtained  by  Calspan  Corporation,  AEDC 
Division,  operating  contractor  of  the  aerospace  flight  dynamics  testing  efforts  at  the  AEDC, 
AFSC,  Arnold  Air  Force  Base,  Tennessee.  The  research  was  performed  in  the  von  K&rm&n 
Gas  Dynamics  Facility  (VKF),  during  the  period  of  October  1,  1986  through  June  26,  1987, 
under  the  AEDC  Project  Number  DB95VW  (Calspan  Number  V32C-CG).  The  manuscript 
was  submitted  for  publication  March  7,  1988  and  was  originally  prepared  in  partial  completion 
of  the  requirements  for  a  Master’s  Degree  at  the  University  of  Tennessee  Space  Institute, 
Tullahoma,  Tennessee. 

The  author  wishes  to  extend  special  thanks  to  Dr.  Bruce  W.  Bomar,  for  his  encouragement, 
suggestions,  and  criticism;  to  Dr.  Julius  Bendat,  for  providing  the  initial  encouragement  to 
apply  the  Hilbert  transform  to  radar  signals;  to  the  Engineering  Applications  Section  of  OAO, 
Incorporated,  for  computer  graphics  support;  and  to  Mr.  Howard  Harris,  for  providing 
countless  technical  details  about  the  Range  G  radar  system  and  data. 
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1.0  INTRODUCTION 
1.1  BACKGROUND  INFORMATION 

A  vehicle  encounters  tremendous  heating  loads  during  atmospheric  reentry.  In  order  to 
survive,  it  must  have  an  effective  heat  shield.  Because  flight  testing  is  very  expensive,  candidate 
heat  shield  mateials  are  often  screened  through  ground  testing. 

One  such  ground  testing  technique  used  at  the  Arnold  Engineering  Development  Center 
(AEDC)  is  to  attach  a  small  specimen  of  candidate  heat  shield  material  onto  a  support  body 
and  launch  the  model  into  a  controlled  atmosphere  at  hypervelocity.  This  procedure  is  often 
used  in  the  von  K&rm&n  Gas  Dynamics  Facility  (VKF)  Hypervelocity  Range/Track  (G). 

A  model  is  accelerated  to  hypervelocity  in  a  short  distance  and,  thus,  encounters  high 
axial  accelerations,  often  in  excess  of  100,000  g.  A  one-dimensional  (1-D)  gas  dynamics  code 
predicts  that  the  problem  of  high-g  loading  is  further  compounded  when  shock  waves  in  the 
launching  gas  reflect  from  the  aft  end  of  the  model  and  cause  acceleration  peaks.  It  has  been 
hypothesized  that  these  peak  accelerations  have  caused  material  sample  failures  during  some 
launches.  (More  recently,  the  high  accelerations  of  launch  have  been  exploited  for  survivability 
testing  of  small  electronic  packages.) 

For  test  planning  purposes,  the  1-D  gas  dynamics  model  is  used  to  predirt  the  velocity 
of  models  launched  in  Range  G.  In  order  to  validate  the  accuracy  of  the  prediction  program, 
it  is  desirable  to  know  the  actual  velocity  and  acceleration  histories  of  models  while  they 
are  still  in  the  gun  barrel.  Thus,  a  Doppler  radar  system  has  been  installed  in  Range  G  in 
order  to  (1)  aid  validation  of  the  launcher  performance  prediction  program,  (2)  assist 
investigations  of  material  sample  failures,  and  (3)  measure  the  accelerations  imparted  to 
electronics  packages. 

The  Doppler  radar  system  produces  a  signal  with  a  frequency  proportional  to  the  model 
velocity.  Previous  attempts  to  calculate  acceleration  from  recorded  radar  signals  have  included 
manually  measuring  the  time  intervals  between  signal  peaks  and  zero  crossings.  Because  the 
peak  interval  measurements  are  labor  intensive,  and  because  the  calculated  results  are  adversely 
affected  by  even  small  errors  in  the  interval  measurements,  these  methods  have  achieved  only 
limited  success.  Thus,  an  alternate  method  of  obtaining  the  model  acceleration  from  the  radar 
signal  was  sought. 

1.2  STATEMENT  OF  THE  PROBLEM  AND  STUDY  OBJECTIVE 

The  acceleration  histories  of  models  launched  in  Range  G  must  be  measured  to  validate 
launcher  performance  prediction  methods  and  to  aid  investigations  of  material  specimen 
failures.  The  objective  of  this  study  was  to  develop  and  verify  a  computer  program  for 
accurately  and  efficiently  calculating  model  accelerations  from  Doppler  radar  data. 
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1.3  APPROACH 

A  literature  search  was  conducted  in  order  to  find  a  better  method  of  calculating  the 
instantaneous  model  accelerations.  The  search  revealed  that  the  Hilbert  transform  was  well- 
suited  for  the  job  because  it  can  be  used  to  obtain  the  instantaneous  frequency  of  a  signal 
and  because  it  can  be  easily  calculated  using  the  fast  Fourier  transform  (FFT).  For  the  case 
of  the  Range  G  Doppler  radar  signal,  the  instantaneous  frequency  is  directly  proportional 
to  the  model  velocity.  It  was  then  determined  that  a  computer  program  could  be  developed 
that  would  obtain  the  model  velocity  from  the  Doppler  radar  signal  via  the  Hilbert  transform 
and  then  numerically  differentiate  the  velocity  to  obtain  the  model  acceleration  history. 

The  Doppler  signal  from  a  typical  model  launch  contains  tens  of  thousands  of  data  points. 
It  was  found  that  the  AEDC  CRAY®  XMP/12  supercomputer  had  available  an  FFT 
subroutine  suitable  for  processing  such  large  arrays  of  data.  The  approach  of  the  study  was 
then  to  implement  the  Hilbert  transform  on  the  CRAY.  A  numerical  differentiation  algorithm 
was  then  added  to  form  the  basis  of  the  radar  data  reduction  program. 

The  program  was  applied  to  simulated  radar  signals.  The  results  were  compared  to  the 
known  theoretical  acceleration  histories.  Random  noise  was  added  to  the  simulated  signals 
to  determine  the  sensitivity  of  the  program  to  imperfect  signals.  The  program  was  then  applied 
to  data  from  several  model  launches,  and  the  results  from  two  radar  systems  were  compared. 

2.0  TEST  FACILITY 

The  Hypervelocity  Range/Track  (G)  shown  in  Fig.  1  is  a  test  unit  that  is  used  for  the 
measurement  of  the  aerodynamic  coefficients  of  models  in  hypervelocity  free  flight  and  for 
captive  testing  of  material  sample  models  along  a  track.  The  facility  can  be  easily  converted 
from  the  free-flight  (range)  configuration  to  the  track  configuration.  The  track  configuration 
allows  the  recovery  of  models  for  posttest  inspection.  The  facility  consists  of  a  launcher, 
a  blast  chamber,  and  the  range/track. 

2,1  LIGHT/GAS  LAUNCHER 

The  Range  G  launcher  (Fig.  2)  is  equipped  with  a  2.5-in.  caliber,  two-stage  powder- 
hydrogen  gun  approximately  170  ft  long.  The  launcher  is  a  fourth-generation  two-stage  gun 
that  has  been  in  use  since  1969  and  has  launched  approximately  3,000  models  at  speeds  up 
to  24,000  ft/sec. 

The  launch  cycle  can  be  summarized  as  follows:  Gun  powder  is  ignited,  which  propels 
a  piston  along  a  pump  tube.  The  piston  compresses  hydrogen  gas  in  the  pump  tube.  When 
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the  hydrogen  pressure  becomes  sufficiently  high,  a  rupture  disk  breaks,  which  then  allows 
the  hydrogen  to  escape  into  the  gun  barrel.  The  high-pressure  hydrogen  then  propels  the 
model  out  the  barrel.  Hydrogen  is  used  because  of  its  high  speed  of  sound.  The  performance 
of  the  launch  cycle  is  controlled  by  varying  the  amount  of  gun  powder,  the  mass  of  the  piston, 
the  initial  pressure  of  the  hydrogen,  the  diaphragm  rupture  pressure,  and  the  mass  of  the  model. 

2.2  INTERIOR  BALLISTIC  DOPPLER  RADAR  SYSTEM 

The  idea  of  using  radar  to  measure  the  interior  ballistics  of  guns  is  not  new.  Doppler 
radar,  or  microwave  interferometry,  was  first  used  to  measure  projectile  motions  in  guns 
as  early  as  the  1950' s  {Refs.  1  and  2).  At  the  AEDC,  a  microwave  interferometer  was  used 
in  Range  C  in  1961  (Refs.  3  and  4). 

The  operation  of  the  Range  G  radar  system  can  be  summarized  as  follows:  A  reflex  klystron 
transmits  microwave  energy  through  an  antenna  that  performs  both  the  transmitting  and 
receiving  functions.  The  microwaves  travel  from  the  barrel  exit  toward  the  model.  The 
microwave  energy  reflects  from  the  front  of  the  model  and  then  returns  to  the  antenna.  The 
reflected  microwaves  are  mixed  with  the  transmitted  microwaves  resulting  in  an  oscillating 
signal  that  is  recorded  on  analog  magnetic  tape.  Equipment  for  operation  at  either  3.6  or 
24  GHz  is  available. 


3.0  THEORETICAL  DEVELOPMENT 


3.1  DOPPLER  RADAR 

When  electromagnetic  waves  are  reflected  from  a  moving  object,  the  reflected  waves 
undergo  a  Doppler  shift  (Ref.  5).  When  sensed  by  an  antenna  that  is  designed  for  sensing 
the  constructive/destructive  interference  of  the  transmitted  and  reflected  waves,  the  antenna 
produces  a  signal  that  undergoes  cyclic  changes  as  the  object  moves.  For  the  case  of  a  model 
in  the  Range  G  launcher,  the  signal  will  go  through  one  complete  cycle  for  every  half¬ 
wavelength  travelled  by  the  model  (Fig.  3).  The  Doppler  frequency,  f(t),  is  described  by  Eq.  (1). 

««)  =  (1) 

where  x(t)  is  model  velocity,  X  is  the  radar  wavelength,  and  t  is  time. 

It  should  be  pointed  out  that  the  wavelength,  X,  is  the  wave  guide  wavelength,  not  the 
free-space  wavelength.  Two  radar  systems  are  considered  in  this  report.  For  the  3.6-GHz 
system,  X  is  0.422  ft;  and  for  the  24-GHz  system,  X  is  0.0424  ft. 
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3.2  PREVIOUS  DATA  REDUCTION  METHODS 

Doppler  radar  has  been  in  use  for  measuring  the  velocity  of  moving  objects  since  the 
early  1940*s  when  it  was  used  to  measure  aircraft  velocity  (Ref.  6).  For  most  applications, 
velocity  is  quasi-steady  (i.e.  low  accelerations),  and  typically  simple  tracking  filters  are  used. 
However,  in  hypervelocity  guns  such  as  the  Range  G  launcher,  the  velocity  and  Doppler 
frequency  change  drastically  in  just  a  few  milliseconds.  In  order  to  get  a  valid  acceleration 
history,  the  Doppler  signal  must  be  analyzed  in  detail  for  an  extremely  large  range  of 
frequencies.  Thus,  past  methods  of  analyzing  Doppler  radar  signals  from  gun  barrels  have 
focused  on  measuring  the  duration  of  each  cycle  of  the  signal. 

The  current  method  used  to  reduce  the  Doppler  radar  data  from  the  Range  G  launcher 
is  based  on  Eq.  (1),  which  can  be  rearranged  as 

x(.)  =  mb  ® 

2 


If  one  defines  the  duration  of  the  ith  half-cycle  of  the  signal  as  AT;,  then  Eq.  (2)  is 
approximated  by 


_X_ 
”  4ATj 


(3) 


A  computer  graphics  terminal  is  used  so  that  an  operator  can  position  a  cursor  over  each 
signal  peak  in  order  to  measure  the  values  of  ATj  for  all  the  half-cycles  of  a  model  launch. 
Hence,  this  technique  is  dubbed  the  “peak-picking”  method.  The  velocity  history  of  Eq. 
(3)  is  then  differentiated  to  obtain  the  model  acceleration. 


3.3  HILBERT  TRANSFORM 


The  Hilbert  transform  was  developed  by  the  German  mathematical)  David  Hilbert 
(1862-1942).  Since  its  development,  the  transform  has  seen  only  limited  application;  and  in 
light  of  the  popularity  of  the  Fourier  transform,  the  Hilbert  transform  can  be  considered 
obscure.  It  does  possess  some  rather  unique  properties  (Refs.  7  and  8)  that  have  led  to  its 
use  in  some  relatively  diverse  fields.  Some  recent  examples  include  image  enhancement  (Ref. 
9),  structural  dynamic  damping  measurements  (Ref.  10),  characterization  of  nonlinearities 
in  structures  (Ref.  11),  and  electronic  communication  (Refs.  12  and  13). 

The  Hilbert  transform  of  an  oscillating  time-wise  function,  y(t),  is  a  second  function  of 
time,  y(t),  where 
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T(t)  =  H[y(t)]  =  {*  du  <4> 

If  we  create  a  complex  analytic  signal  (Refs.  7  and  8), 

z(t)  =  y(t)  +  jyft)  ,  j  =  V3!  (5) 

we  can  obtain  some  very  useful  properties  as  shown  in  the  following  text. 

Any  complex  variable  can  be  written  in  exponential  form,  so  z(t)  can  be  written  as 

z(t)  =  A(t)exp[j</>(t)]  (6) 

where 

A(t)  =  Wt)  +  y!(t)]1/!  (7) 

and 

*(t)  =  Tan->[|$]  (8) 

Thus,  z(t)  is  a  helix  in  complex  space  with  radius  A(t)  and  phase  angle  <$(t)  as  shown  in  Fig. 

4.  A(t)  is  also  the  instantaneous  amplitude  of  y(t)  and  of  y(t).  The  instantaneous  frequency 
is  given  by 

«.)  =  f  l9) 

from  which  we  can  calculate  velocity  using  Eq.  (2).  Thus,  if  the  Hilbert  transform  of  the 
Range  G  Doppler  radar  signal  can  be  obtained,  the  model  position,  velocity,  and  acceleration 
histories  can  be  easily  calculated. 

A  method  for  obtaining  the  Hilbert  transform  of  a  signal  via  the  Fourier  transform  and 
the  analytic  signal  is  presented  in  Ref.  8  and  is  shown  here. 

The  Fourier  transform  is  defined  as 

Y(f)  =  F[y(t)]  =  J"  y(t)exp  (— j2irft)  dt  (10) 

—  as 

The  convolution  property  of  the  Fourier  transform  gives 
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F[yi(t)  *  y2(t)]  =  F[y i (t)] F [y2(t>]  (11) 

=  Yi(0y2(0 

where  *  denotes  convolution,  which  is  defined  as 

yi(0  *  y2(0  =  pyiO^G  -  u)  du  (12) 


Since  the  definition  of  the  Hilbert  transform  [Eq.  (4)]  is  equivalent  to  the  convolution 
of  y(t)  and  i/rt;  i.e., 

7(0  =  y(0  •  s  03) 

We  can  use  Eq.  (11)  to  write 


where 


F[  y(t)]  =  F^y(t)  *  J  =  FtyftJjF^) 

=  Y(f)[  -j  •  sgn(f)] 


sgn(f)  = 


I 


-1,  f  <  0 
1,  f  >  0 


(14) 


(15) 


The  Fourier  transform  of  the  analytic  signal  [Eq.  (5)]  is 

Z(f)  =  F[z(t)J  =  F[y(t)]  +  jF[  y(t)]  (16) 

Upon  substituting  Eq.  (14),  we  obtain 


where 


Z(f)  =  Y(f)  +  j Y(f)[  -  j  .  sgn(f)] 

=  Y(f)[l  +  sgn(f)]  =  B(f)Y(f) 


B(f) 


0,  f  <  0 
2,  f  >  0 


(17) 


(18) 
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From  inspection  of  Eqs.  (17)  and  (18),  we  can  see  that,  given  the  Fourier  transform  of 
y(t),  we  very  easily  obtain  the  Fourier  transform  of  z(t)  by  setting  Y(f)  to  zero  for  negative 
frequencies  and  by  multiplying  Y(f)  by  2  for  positive  frequencies.  The  analytic  signal  can 
then  be  obtained  via  the  inverse  Fourier  transform. 

3.4  COMPARISON  OF  TWO  METHODS 

The  primary  advantage  of  the  analytic  signal  method  of  calculating  velocity  over  the  peak¬ 
picking  method  is  illustrated  in  Fig.  5.  The  analytic  signal  method  provides  a  continuous 
estimate  of  the  signal  phase  angle,  whereas  the  peak-picking  method  provides  a  stair-step 
phase  angle  history.  Thus,  the  differentiation  needed  to  obtain  the  model  velocity  can  be 
more  readily  performed  using  the  continuous  curve  of  the  analytic  signal  method  rather  than 
the  disjointed  curve  of  the  peak -picking  method. 

3.5  SMOOTHING  AND  DIFFERENTIATION 

Numerical  differentiation  of  a  noisy  function  typically  magnifies  the  effect  of  the  noise. 
If  one  can  tolerate  a  loss-of-time  resolution,  numerical  smoothing  can  be  applied  either  before 
or  after  the  differentiation  to  reduce  the  effect  of  the  magnified  noise.  The  Hanning  filter 
is  a  well-known  smoothing  technique  and,  as  shown  below,  can  be  combined  with  the  numerical 
differentiation  into  a  single  operation. 

The  normalized  Hanning  window,  h(t),  is  defined  as 


(19) 


where  T  is  the  time  span  of  the  window.  The  low-pass  nature  of  Hanning  smoothing  is 
illustrated  in  Fig.  6. 

If  x(t)  is  a  function  to  be  differentiated,  and  x(t)  is  the  result  of  the  differentiation,  x(t) 
can  be  smoothed  by  convolution  with  the  Hanning  window  yielding 
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x'(t)  =  h(t)  *  x(t)  =  J“h(u)x(t  -  u)  du 


(20) 


Integration  by  parts  gives 


x'(0  -  -  h(u)x(t  -  u)|  +|  h(u)x(t  -  u) 

II  =  —  oo  —  oo 


du 


(21) 


Substitution  of  Eq.  (19)  gives 


T 

u  =  — 

2  2* 


*'<o  =  °+  I  t  S sin (¥)x(t - U>  dU 


(22) 


This  is  equivalent  to  convolution  with  a  differentiated  Hanning  window;  i.e., 

x'(t)  =  h(t)  *  x(t)  =  h(t)  *  x(t)  (23) 

Thus,  we  can  obtain  a  smoothed,  differentiated  version  of  x(t)  by  a  single  convolution  with 
h(t),  where 


h, 


h(t)  = 


sin  2nrt, 
T2  T" 


t  <  ~T 
2 

^T<  t  <  T 
2  2 


0, 


t  > 


T 

2 


(24) 


This  technique  is  simple  to  implement  numerically  and  is  used  to  calculate  4>(t)  and 
acceleration  as  discussed  in  Section  4.5. 


4.0  PROGRAM  DESCRIPTION 

A  Fortran  program  called  TXHARS  (Trajectory  Extraction  via  Hilbert  Analysis  of  a  Radar 
Signal)  has  been  written  using  the  techniques  described  in  Section  3.  An  interactive  program 
that  runs  on  the  AEDC  AMDAHL  computer  queries  the  user  for  the  necessary  parameters. 
The  query  program  then  automatically  submits  the  TXHARS  program  to  the  CRAY  for 
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execution.  After  TXHARS  executes,  the  results  are  sent  back  to  the  AMDAHL  where  the 
user  can  view  the  results  using  TEKPLOT,  an  interactive  graphics  program  available  at  the 
AEDC. 

4.1  PROGRAM  SUMMARY 

Detailed  descriptions  of  the  various  algorithms  used  in  the  TXHARS  program  are  presented 
in  Sections  4.2  through  4.7.  However,  the  basic  operation  of  the  program  can  be  summarized 
as  follows: 

1.  A  helix  in  complex  space  is  created  via  the  Hilbert  transform. 

2.  Band-pass  filtering  is  performed  in  order  to  remove  noise. 

3.  The  phase  angle,  <f>( t),  of  the  helix  is  simultaneously  smoothed  and  differentiated 
to  obtain  the  model  velocity,  x(t). 

4.  A  correction  is  applied  to  the  computed  velocity  during  periods  of  insufficient 
signal  amplitude. 

5.  The  corrected  velocity  is  integrated  to  obtain  the  model  position  history,  x(t). 

6.  The  corrected  velocity  is  simultaneously  smoothed  and  differentiated  to  obtain 
the  model  acceleration,  x(t). 

4.2  CALCULATING  THE  HILBERT  TRANSFORM 

For  the  real-world  case  of  finite  signal  lengths  and  discrete  data  points,  the  procedure 
presented  in  Section  3.3  for  obtaining  the  Hilbert  transform  can  be  approximated  using  the 
discrete  Fourier  transform  (DFT)  as  calculated  by  the  fast  Fourier  transform  (FFT)  algorithm. 
The  negative  frequencies  of  Eqs.  (15)  and  (18)  are  manifested  as  frequencies  between  the 
half-sample  frequency  and  the  sample  frequency.  The  resulting  procedure  for  obtaining  the 
analytic  signal  is  as  follows: 

1.  Compute  the  FFT  of  y(t)  to  obtain  Y(f). 

2.  Multiply  Y(f)  by  2  for  frequencies  betwen  zero  and  the  half-sample  frequency. 
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3.  Set  Y(f)  to  zero  for  frequencies  between  the  half-sample  frequency  and  the  sample 
frequency.  [Steps  2  and  3  yield  Z(f)-1 

4.  Compute  the  inverse  FFT  of  Z(f)  to  obtain  z(t). 

The  TXHARS  program  uses  the  IMSL  subroutine  FFTCC  (Ref.  14)  to  compute  the  FFT 
of  the  radar  signal  and  the  inverse  FFT  of  the  analytic  signal  to  obtain  the  Hilbert  transform 
of  the  radar  signal.  Since  the  subroutine  does  not  have  a  constraint  on  the  number  of  data 
points,  the  very  large  array  defining  the  radar  signal  can  be  processed  with  only  one  call 
for  the  FFT  and  a  second  call  for  the  inverse  FFT.  This  is  simpler  and  more  efficient  than 
segmenting  the  signal  into  several  short  pieces  and  making  several  calls  to  the  FFT  subroutine. 
To  date,  radar  signals  containing  up  to  40,000  data  points  have  been  processed. 

4,3  SIGNAL  FILTERING 

Since  the  signals  from  the  radar  system  are  typically  noisy,  provisions  for  filtering  the 
radar  signal  are  included  in  TXHARS.  Since  signal  frequency  is  proportional  to  model  velocity, 
any  signal  content  at  frequencies  higher  than  that  of  the  maximum  velocity  achieved  by  the 
model  can  be  assumed  to  be  noise.  Also,  very-low-frequency  data  (such  as  signal  drift)  can 
also  be  ignored  since  the  model  achieves  a  finite  velocity  in  a  very  short  time  when  hit  by 
the  first  hydrogen  pulse.  Thus,  the  program  user  is  allowed  to  select  a  velcoity  range  of  interest 
that  is  defined  by  VM1N  and  VMAX.  Signal  frequencies  corresponding  to  velocities  outside 
the  range  of  interest  are  eliminated  by  band-pass  filtering. 

The  band-pass  filtering  is  performed  in  the  frequency  domain  by  zeroing  out  Z{f)  at 
frequencies  outside  the  band-pass  limits  f,  and  f2.  This  results  in  a  simple  “box  car”  filter 
as  illustrated  in  Fig.  7.  The  band-pass  limits  are  calculated  by 

f,  -  0.9  f2  „  U  (UMAS)  (25) 

The  factors  outside  of  the  parentheses  serve  to  slightly  broaden  the  pass  band. 

One  shortcoming  of  the  box  car  filter  is  the  phenomenon  known  as  “ripple”  (Ref.  15), 
which  is  illustrated  as  the  scalloped  curve  in  Fig.  7.  Ripple  is  the  result  of  approximating 
a  continuous  spectrum  with  discrete  points.  This  causes  nonzero  filter  components  to  exist 
within  the  stop  bands  and  nonunity  components  to  exist  in  the  pass  band.  A  more  advanced 
filter  that  has  less  ripple  than  the  box  car  could  have  been  used,  but  the  box  car  is  very  easy 
to  implement  and  has  provided  acceptable  results  to  date. 
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After  the  box  car  filter  is  applied  to  Z(f),  we  are  left  with  Z '  (f).  The  inverse  FFT  of  Z  ’  (f) 
then  gives 

z'(t)  -  y'(t)  +  j  y'  (t)  (26) 

where  y'(t)  is  the  filtered  radar  signal,  and  y(t)  is  the  Hilbert  transform  of  y'(t). 

4.4  UNWRAPPING  THE  PHASE  ANGLE 

The  amplitude,  A{t),  and  phase  angle,  of  the  filtered  signal  are  calculated  using  Eqs. 
(7)  and  (8).  However,  since  the  inverse  tangent  function  available  in  Fortran  only  provides 
results  ranging  from  -  ?r  to  +r,  the  <£(t)  values  must  be  “unwrapped”  as  illustrated  in  Fig. 
8.  The  program  has  an  unwrapping  algorithm  that  appropriately  adds  or  subtracts  2s  from 
the  phase  angle  whenever  two  consecutive  points  are  in  the  second  and  third  quadrants. 

4.5  HANNING  SMOOTHING 

TXHARS  uses  the  Hanning  window  methods  presented  in  Section  3.5  to  smooth  the 
velocity,  acceleration,  and  signal  amplitude  histories.  A  window  that  covers  a  large  time  span 
is  effective  at  removing  unwanted  noise  from  the  results,  but  it  also  removes  desired  detail. 
Conversely,  a  window  that  covers  a  small  time  span  provides  better  detail,  but  can  also  pass 
too  much  noise.  Thus,  the  program  was  designed  so  that  the  user  can  adjust  the  time  span 
of  the  windows  in  order  to  reach  an  acceptable  compromise  between  noise  and  detail. 

Since  each  of  the  data  arrays  to  be  smoothed  is  of  finite  length,  the  Hanning  window 
cannot  function  properly  at  the  array  ends.  This  is  because  there  are  no  data  points  available 
for  convolution  with  the  Hanning  window  outside  of  an  array.  Therefore,  for  these  periods 
the  program  uses  the  convolution  value  obtained  at  the  most  extreme  point  for  which  full 
overlap  occurs.  This  results  in  a  small  flat  at  each  end  of  a  smoothed  array.  This  effect  is 
illustrated  in  Fig.  9. 


Figure  10  shows  the  effect  of  Hanning  windows  of  various  time  spans  on  a  rectangular 
pulse.  The  figure  is  presented  in  dimensionless  form.  The  pulse  begins  at  time  zero,  immediately 
obtains  a  value  of  one,  and  then  drops  back  to  zero  at  time  one.  The  R  =  0.5  curve  shows 
that  a  very  good  approximation  of  a  pulse  is  obtained  when  a  rectangular  pulse  is  smoothed 
by  a  Hanning  window  with  a  time  width  one-half  the  pulse  width.  A  very  poor  approximation 
is  obtained  when  the  Hanning  width  is  10  times  the  pulse  width  as  shown  by  the  R  =  10  curve. 

As  an  example,  suppose  one  wanted  to  know  the  peak  g  loading  encountered  during 
acceleration  pulses  lasting  more  than  0. 1  msec.  The  R  =  1  curve  of  Fig.  10  shows  these  peak 
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values  could  be  obtained  with  a  Hanning  window  as  large  as  0. 1  msec.  Note  that  the  detailed 
shape  of  the  acceleration  curve  would  be  incorrect. 

It  should  also  be  pointed  out  that  the  response  curves  of  Fig.  10  apply  to  the  difference 
between  the  pulse  height  and  the  nominal  value,  not  the  absolute  height.  As  another  example, 
suppose  the  nominal  acceleration  of  the  model  was  100,000  g,  then  suddenly  jumped  to  150,000 
g  for  0.1  msec.  The  R  =  2  curve  shows  that  a  Hanning  window  of  a  0.2-msec  span  would 
obtain  the  100,000-g  nominal  acceleration  plus  about  84  percent  of  the  50,000-g  difference. 
Thus,  the  TXHARS  output  for  the  peak  acceleration  would  be  142,000  g.  This  would  be 
considered  a  reasonable  engineering  estimate  of  peak  acceleration  since  it  respresents  an  error 
of  only  5  percent. 

Hanning  windows  with  time  widths  greater  than  the  pulse  width  of  interest  will  always 
underestimate  the  peak  acceleration  value,  whereas  signal  noise  typically  causes  the  program 
to  overestimate  the  peak  values.  Knowledge  of  the  Hanning  window  widths  and  the  signal 
quality  are,  therefore,  useful  when  interpreting  the  program  output. 

4.6  VELOCITY 

An  initial  estimate  of  the  model  velocity,  x{t),  is  found  by  first  applying  Eq.  (24)  to 
simultaneously  smooth  and  differentiate  </>(t)  to  get  $(t);  then  Eqs.  (9)  and  (2)  are  used  to 
get  x(t).  The  velocity,  x(t),  must  then  be  corrected  because  the  amplitude  of  Range  G  radar 
signals  often  drops  to  zero.  Typical  drops  in  amplitude  are  shown  in  Fig.  1 1  and  are  believed 
to  be  caused  by  several  phenomena  that  include  antenna  vibration  and  extraneous  radar 
reflections.  Since  these  amplitude  drops  can  significantly  degrade  the  program  results, 
provisions  are  included  in  TXHARS  to  interpolate  the  velocity  across  these  time  periods. 
This  effectively  “skips”  the  “bad”  portions  of  the  signal.  The  skipping  algorithm  uses  a 
floating  criterion  for  identifying  bad  periods. 

The  floating  criterion  used  in  the  ratio  of  a  lightly  smoothed  trace  of  the  signal  amplitude, 
A(t),  with  a  heavily  smoothed  trace.  Equation  (19)  is  convolved  with  the  data  to  perform 
the  smoothing.  When  the  lightly  smoothed  curve  is  less  than  a  user-specified  fraction  of  the 
heavily  smoothed  curve,  velocity  interpolation  is  performed.  After  interpolation  is  performed, 
Eq.  (19)  is  convolved  with  x(t)  in  order  to  eliminate  kinks  in  the  curve.  The  algorithm  is 
illustrated  in  Fig.  12.  The  effect  of  turning  off  the  skipping  algorithm  is  shown  in  Fig.  13. 
Note  the  large  spikes  on  the  “Without  Velocity  Correction”  velocity  curve.  These  spikes 
correspond  to  signal  amplitude  dropouts. 
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4.7  POSITION  AND  ACCELERATION 

After  the  skipping  algorithm  processes  the  velocity  history,  x(t),  the  position,  x(t),  is 
calculated  by  numerically  integrating  velocity.  A  simple  rectangular  integration  method  is 
used.  Acceleration,  x(t),  is  calculated  by  convolving  the  differentiated  Hanning  window  of 
Eq.  (24)  with  x(t). 

4.8  QUERY  PROGRAM 

Because  several  parameters  can  be  adjusted  by  the  program  user,  an  interactive  query 
program  has  been  added  to  TXHARS  by  the  Engineering  Applications  Section  of  OAO, 
Incorporated,  the  computer  services  subcontractor  at  the  AEDC.  A  typical  session  is  presented 
in  Appendix  A  along  with  definitions  for  the  input  variables.  The  query  program  was  designed 
so  that  the  user  can  easily  track  the  variables  used  to  produce  the  results  of  numerous  launches 
and  analyses.  The  query  program  automatically  produces  a  hardcopy  of  the  terminal  screen 
at  the  end  of  each  session.  The  session  date  and  time  are  included  on  the  hardcopy  and  stored 
for  subsequent  labeling  of  the  plotted  results. 

5.0  RESULTS 


5.1  SIMULATED  SIGNALS 

Several  simulated  signals  have  been  generated  in  order  to  validate  the  TXHARS  computer 
program.  The  acceleration  histories  produced  by  TXHARS  were  compared  to  the  theoretical 
accelerations  corresponding  to  each  of  the  simulated  signals.  Figure  14  shows  the  theoretical 
model  trajectory  for  one  such  simulation.  Figure  15  shows  the  velocity  and  acceleration  histories 
calculated  by  TXHARS  for  a  simulated  signal  with  a  Doppler  wavelength  of  the  3.6-GHz 
radar  system,  and  Fig.  16  shows  the  trajectory  calculated  for  a  simulated  24-GHz  signal. 
The  close  agreement  of  the  calculated  trajectories  in  Figs.  15  and  16  with  the  theoretical 
trajectory  of  Fig.  14  validates  the  basic  operation  of  the  program. 

It  can  be  seen  that  the  acceleration  curves  of  the  simulated  signals  are  oscillatory  in  nature. 
The  period  of  the  oscillation  was  selected  to  correspond  to  a  duration  of  acceleration  peak 
of  0.08  msec.  This  duration  is  typical  of  the  acceleration  pulses  predicted  by  the  launcher 
performance  program. 

It  can  be  seen  from  Figs.  14  and  15  that  both  the  3.6-  and  24-GHz  signals  could  be  used 
to  resolve  these  peaks  reasonably  well.  Acceleration  peaks  of  shorter  duration  (i.e.  higher 
frequency  oscillations  of  the  theoretical  acceleration)  could  not  be  resolved  nearly  as  well, 
with  the  results  tending  to  be  worse  at  the  lower  theoretical  model  velocities,  especially  for 
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the  3.6-GHz  runs.  Thus,  0.08  msec  can  be  used  as  a  rule  of  thumb  as  the  smallest  acceleration 
peak  resolvable  by  TXHARS  with  perfect  signal  quality. 

In  order  to  investigate  the  effect  of  signal  noise  on  the  program  performance,  various 
amounts  of  band-limited  white  noise  were  added  to  a  simulated  signal.  The  random  noise 
was  generated  using  the  IMSL  subroutine  GGNML  (Ref.  14).  The  theoretical  trajectory  chosen 
for  the  noise  study  was  similar  to  the  one  in  Fig.  14,  except  the  acceleration  oscillation  period 
was  increased  to  represent  an  acceleration  peak  duration  of  0.1  msec.  Figure  17  shows  the 
TXHARS  results  for  a  3.6-GHz  signal  with  increasing  noise  content,  and  Fig.  18  shows  the 
24-GHz  results. 

Note  that  the  velocity  histories  are  almost  unaffected  by  noise,  whereas  the  acceleration 
histories  deviate  significantly  with  increased  noise.  The  results  can  be  summarized  by  stating 
that  an  acceleration  peak  of  0. 1  msec  can  be  reliably  resolved  with  a  signal-to-noise  ratio 
as  low  as  30  db  for  the  3.6-GHz  signals  and  20  db  for  24-GHz  signals.  Acceleration  peaks 
of  longer  duration  can  be  resolved  with  even  lower  signal-to-noise  ratios.  These  results 
reemphasize  that  one  needs  to  consider  signal  quality  when  evaluating  the  program  results 
for  a  specific  launch. 

5.2  RANGE  G  RADAR  SIGNAL 

The  Doppler  radar  signals  from  several  model  launches  have  been  analyzed  with  the 
TXHARS  program,  and  typical  results  are  presented  in  Figs.  19  and  20.  These  analyses  were 
made  using  the  following  parameter  values,  which  are  defined  in  Appendix  A: 

CUT  =  0.7 
TFILT1  =  0.20  msec 
TFILT2  =  0.05  msec 
TFILT3  =  0.30  msec 
TFILT4  =  0.40  msec 

Since  the  true  model  trajectories  are  not  known,  absolute  statements  about  the  program 
accuracy  cannot  be  made  based  on  these  results.  However,  some  interesting  trends  can  be 
cited  as  noted  in  the  following. 
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As  mentioned  in  Section  4.6,  the  amplitude  of  the  radar  signal  occasionally  drops  to  near 
zero.  Theoretically,  this  should  not  pose  a  problem  since  the  phase  angle  of  the  analytic  signal 
should  continue  to  increase  during  these  periods.  However,  Fig.  21  shows  that  for  actual 
data,  the  helix  of  the  analytic  signal  typically  is  not  centered  around  the  time  axis  during 
the  signal  amplitude  drop.  This  offset  can  be  attributed  to  low-frequency  noise  in  the  signal. 
The  offset  causes  the  phase  angle  to  oscillate,  instead  of  increase,  with  time.  The  velocity 
correction  algorithm  described  in  Section  4.6  has  proven  to  be  successful  in  eliminating  the 
adverse  effects  of  the  amplitude  drop,  but  an  improved  filtering  method  might  provide  better 
results.  Such  a  method  would  have  to  eliminate  the  helix  offset  of  the  amplitude  drops  without 
filtering  out  the  low-frequency  Doppler  signal  that  occurs  during  the  early  period  of  a  model 
launch. 

It  is  believed  that  one  measure  of  signal  quality  is  the  size  of  the  Hanning  windows  used 
for  smoothing  the  data.  Typically,  when  a  launch  signal  is  analyzed  for  the  first  time,  very 
small  Hanning  windows  are  selected,  and  usually  very  noisy  and  unrealistic  results  are 
produced.  The  program  is  then  reexecuted  with  increasing  window  widths  until  the  program 
results  are  judged  to  be  reasonable.  (Usually  no  more  than  three  or  four  runs  are  needed 
to  get  acceptable  results.) 

Signals  from  the  3.6-GHz  system  often  require  smaller  windows  and,  thus,  are  felt  to 
produce  better  results  than  the  24-GHz  system.  This  would  seem  to  be  contrary  to  the  results 
of  the  simulated  signal  study  of  Section  S.l  that  indicate  that  the  higher  24-GHz  system  should 
provide  better  results.  However,  experience  has  shown  that  the  24-GHz  system  produces  much 
noisier  signals  than  the  3.6-GHz  system  and  that  the  amplitude  of  the  signal  from  the  3.6-GHz 
system  is  less  prone  to  drop  to  zero.  In  practice,  therefore,  the  inferior  signal  quality  of  the 
24-GHz  system  usually  negates  any  inherent  advantage  it  has  over  the  3.6-GHz  system. 

5.3  PERFORMANCE 

Signals  of  various  length  have  been  analyzed  with  TXHARS,  the  largest  containing  40,000 
data  points.  Figure  22  shows  the  adverse  effect  of  signal  length  on  run  time.  However,  the 
maximum  central-processor  time  of  20  sec  is  considered  to  be  a  small  job  for  the  AEDC 
CRAY.  The  ellapsed  time  from  job  submittal  to  plotted  results  was  typically  13  min  and 
was  largely  unaffected  by  signal  length. 

Since  the  computation  time  of  the  FFT  is  proportional  to  N  •  log(N)  (Ref.  14),  where 
N  is  the  signal  length,  we  can  conclude  from  the  linearity  of  Fig.  22  that  the  run  time  of 
TXHARS  is  not  significantly  influenced  by  the  FFT  calculation.  This  suggests  that  data  input 
and  output  are  the  major  contributors  to  run  time.  Improved  file  handling  methods  might 
serve  to  make  the  program  even  more  efficient. 
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5.4  FUTURE  WORK 

An  effort  is  currently  underway  at  the  AEDC  to  acquire  data  simultaneously  with  both 
the  3.6-  and  24-GHz  systems.  Also,  these  data  will  be  analyzed  using  both  the  TXHARS 
program  and  the  peak-picking  methods.  The  results  of  this  study  should  help  reveal  the 
strengths  and  weaknesses  of  the  different  systems  and  methods.  New  methods  will  also  be 
introduced  in  an  attempt  to  increase  radar  signal  quality. 

A  new  data  acquisition  system  is  currently  being  implemented  in  the  Range  G  facility. 
This  data  system  will  have  the  capability  to  directly  digitize  the  signal  from  the  interior  ballistics 
radar  system  as  a  model  is  launched.  Thus,  the  need  for  magnetic  tape  will  be  eliminated. 
Since  the  internal  memory  capacity  of  the  data  acquisition  system  computer  is  several  times 
larger  than  needed  to  hold  the  radar  signal  history  from  a  typical  launch,  it  may  be  feasible 
to  run  the  TXHARS  program  on  the  new  computer.  This  would  be  desirable  since  it  would 
eliminate  the  data  file  transfers  presently  required  to  move  Range  G  radar  signals  to  the  CRAY. 
TXHARS  could  then  provide  results  more  promptly  if  it  operated  on  the  new  data  acquisition 
system  computer. 

In  order  to  implement  TXHARS  on  the  new  data  acquisition  computer,  two  obstacles 
must  be  overcome.  First,  an  alternate  method  of  calculating  the  Hilbert  transform  would 
have  to  be  used  since  a  suitable  FFT  routine  is  not  available  on  this  computer  [One  alternative 
is  a  finite-impulse-response  digital  filter  (Ref.  15).]  Second,  TXHARS  would  have  to  be 
modified  to  operate  with  the  plotting  software  available  on  the  new  computer,  and  it  is  not 
certain  if  the  plotting  software  is  suitable.  Thus,  further  investigation  is  needed  to  determine 
if  it  is  warranted  to  convert  TXHARS  for  use  on  the  new  data  acquisition  system  computer. 

6.0  CONCLUSIONS  AND  RECOMMENDATIONS 

The  objective  of  this  study  was  met  by  developing  a  data  reduction  program  for  the  Range 
G  interior  ballistics  radar  system.  The  Hilbert  transform  has  been  shown  to  be  well-suited 
for  extracting  model  ballistic  information  from  Doppler  radar  signals.  The  program  has  been 
shown  to  be  accurate  and  efficient.  Signals  from  the  3.6-GHz  system  have  proved  to  provide 
better  results  than  those  from  the  24-GHz  system. 

It  is  recommended  that  the  use  of  more  advanced  filtering  methods  be  pursued  in  order 
to  improve  the  ability  of  the  program  to  deal  with  noisy  signals.  Data  acquired  simultaneously 
with  both  the  3.6-  and  24-GHz  systems  should  be  evaluated,  efforts  to  improve  the  quality 
of  the  radar  signals  from  future  tests  in  Range  G  should  be  continued,  and  the  feasibility 
of  implementing  TXHARS  on  the  new  data  acquisition  system  computer  should  be  studied 
further. 
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Figure  2.  Range  G  launcher. 
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Figure  5.  Analytic  signal  and  peak-picking  methods. 
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Figure  8.  Phase  angle  unwrapping. 
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Figure  9.  Typical  end  effect  of  smoothing  routine. 
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Figure  10.  Effect  of  Hanning  smoothing  on  square  pulse. 
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Figure  12.  Velocity  correction  algorithm. 
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Figure  13.  Effectiveness  of  velocity  correction  algorithm 
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Figure  IS.  TXHARS  results  with  perfect  3.6-GHz  signal. 
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Figure  16.  TXHARS  results  with  perfect  24-GHz  signal. 
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Time,  t,  msec  : 

a.  40-db  signal-to-noise  ratio 
Figure  17.  Effect  of  noise  on  simulated  3.6-GHz  results. 
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b.  30-db  slgnal-to-noise  ratio 
Figure  17.  Concluded. 
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a.  30-db  signal-to-noise  ratio 
Figure  18.  Effect  of  noise  on  simulated  24-GHz  results. 
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c.  15-db  signal-to-noise  ratio 
Figure  18.  Concluded. 
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APPENDIX  A 

QUERY  PROGRAM  PARAMETERS  AND  OPERATION 

The  program  user  will  be  prompted  for  the  following  variables: 

N  The  number  of  data  points  to  be  included  in  the  analysis.  Typically,  numerous 

data  points  before  and  after  the  actual  launch  are  recorded  and  stored  on  the 
computer.  The  user  should  select  only  those  points  that  include  the  model  launch. 

LAMBDA  The  effective  radar  wavelength  (in  feet)  in  the  barrel. 

CUT  The  ratio  of  the  lightly  smoothed  to  heavily  smoothed  amplitude  histories  used 
by  the  skipping  algorithm  for  establishing  a  bad  time  period.  Velocity 
interpolation  across  bad  periods  is  performed.  A  value  of  zero  for  CUT  turns 
off  the  skipping  algorthim. 

TFILT1  The  time  span  (in  milliseconds)  of  the  differentiated  Hanning  window  used  for 
computing  the  uncorrected  velocity  history. 

TFILT2  The  time  span  (in  milliseconds)  of  the  Hanning  window  used  for  the  lightly 
smoothed  amplitude  trace.  The  skipping  algorithm  multiplies  this  number  by 
7  to  obtain  the  size  of  the  Hanning  window  for  computing  the  heavily  smoothed 
amplitude  trace. 

TFILT3  The  time  span  (in  milliseconds)  of  the  Hanning  window  used  for  smoothing  the 
“corrected”  velocity  history  produced  by  the  skipping  algorithm. 

i 

TFILT4  The  time  span  (in  milliseconds)  of  the  differentiated  Hanning  window  used  for 
computing  the  acceleration  history. 

VMIN  The  minimum  velocity  (in  feet/second)  of  interest.  Used  for  establishing  the  lower 
pass-band  frequency  for  filtering  of  the  radar  signal. 

VMAX  The  maximum  velocity  (in  feet/second)  of  interest.  Used  for  establishing  the 
upper  pass-band  frequency  for  filtering  of  the  radar  signal. 

TSTART  Time  (in  milliseconds)  of  the  first  point  of  the  stored  data  to  be  included  in  the 
analysis. 
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Sample  1.  Typical  Query  Session 

*nt.  it.  in?  am  p.h. 


ijtcj  nunieo  or  points  in  siq.  (default  ■  40000 > 

J  *vvv 

N  30000 


ENTER  RADAR  WAVELENGTH  IN  FttT.  (DEFAULT  ■  0.04) 
0.04334 

LAMM  *  0.043 


ENTER  FRACTION  FOR  FUMTtNO  CRITCRON  FOR  UCLOCITV  CORRSCTXON.  (DEFAULT  •  0.7) 

CUT  •  0.700 

ENTER  SIZE  OF  NANNING  WINDOW  (HXLLXSEC)  FOR  1ST  UCLKITV  CALCULATION  I  DEFAULT  -  0.03) 
TFILT1  •  0.100 

ENTER  SIZE  OF  NANNING  WINDOW  tNILLISCC)  FOR  AMPLITUDE  SNOO THING  (DEFAULT  •  0.03) 
TFILT8  -  0.100 

ENTER  SIZE  OF  HANNINO  WINDOW  (ML  LI  SEC)  FOR  UEL0C1TV  SNOO  THING  (DEFAULT  ■  0.03) 

TFILT3  -  0.10# 

ENTER  SIZE  OF  HANNINO  WINDOW  (NXLLIMC)  FOR  ACCELERATION  CALCULATION  (DEFAULT  •  0.03) 
TFILT4  -  0.100 

ENTER  MINIMS!  VELOCITY  IN  FT/SCC  (DEFRUIT  •  100.) 

WHIN  •  800.000 


ENTER  NAM  RUN  VELOCITY  JM  FTsOEC  (DEFAULT  •  10000.  > 
MW  -  10000.000 

ENTER  OTART  TIRE  IN  RXLLI080  (DEFAULT  -  7S.0) 

BO. 

T ST ART  .  80.000 
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NOMENCLATURE 

A  Amplitude  of  analytic  signal 

F[  ]  Fourier  transform 

f  Frequency,  Hz 

fi  Lower  band-pass  limit  of  box  car  filter,  Hz 

fj  Upper  band-pass  limit  of  box  car  filter,  Hz 

H[  1,  ~  Hilbert  transform 

h  Hanning  window,  sec-1 

j 

N  Total  number  of  points  in  signal 

R  Ratio  of  acceleration  pulse  duration  to  Hanning  window  time-span 

T  Time-span  of  Hanning  window,  sec 

ATj  Duration  of  ith  half-cycle  of  radar  signal,  sec 

t  Time,  sec 

(t)  Instantaneous  function  of  time 

u  Dummy  integrand 

x  Model  position,  ft 

x;  Average  model  velocity  during  ith  half-cycle  of  radar  signal,  ft/sec 

Y(f)  Fourier  transform  of  radar  signal 

y  Measured  radar  signal 
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yi  Arbitrary  function  Number  1 

y2  Arbitrary  function  Number  2 

Z(f)  Fourier  transform  of  analytic  signal 

z  Analytic  signal 

'  Low-pass  filter 

*  Convolution 

First  derivation  with  respect  to  time 
Second  derivation  with  respect  to  time 
X  Radar  wavelength  in  gun  barrel,  ft 

<f>  Phase  angle  of  analytic  signal,  radians 
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